Inverse Linear Regression

Inverse Linear Regression

Introduction

The goal:  fit bi-variate data (x,y) to the curve:

y = 1 / (a*x + b)

We will be able to use the linear regression model with the following transformation:

1/y = a* x + b

y' = 1/y
x' =  x
a' = a
b' = b

y' = a' * x' + b'

Let's illustrate this with an example.

Example

Fit the curve y = 1 / (a*x + b) to the data:

(-2, -1.43)
(-1, 0.4)
(0, 0.18)
(1, 0.11)
(2, 0.08)
(4, 0.05)

Transform the data to (x', y'): x' = x,  y' = 1/y

(-2, -1/1.43 = 0.6693006993)
(-1, 1/0.4 = 2.5)
(0, 1/0.18 = 5.555555556)
(1, 1/0.11 = 9.090909091)
(2, 1/0.08 = 12.5)
(4, 1/0.05 = 20)

Regression analysis with the transformed data:

Slope (a,m) = 3.443917194
Intercept (b) = 5.861915862
r^2 = 0.998386787

1/y = 3.443917194 * x + 5.861915862

y = 1 / (3.443917194 * x + 5.861915862)

Note:  My next blog entry will be posted on September 2, 2019, 12:00 AM Pacific Daylight Time on Labor Day. 

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

HP Prime and TI Nspire CX CAS: Solving Integral Equations

HP Prime and TI Nspire CX CAS:  Solving Integral Equations

Introduction

The program INTEGRALSOLVE solve the following equation:
(Format:  ∫( integrand dvar, lower, upper)

∫( f(t) dt, 0, x) = a

∫( f(t) dt, 0, x) - a = 0

It is assumed that x>0. 

We can use the Second Theorem of Calculus which takes the derivative of the integral:

d/dx  ∫( f(t) dt, a, x) = f(x)

We don't have to worry about lower limit a at all for the theorem to work. 

∫( f(t) dt, 0, x) - a

Take the derivative with respect to x on both sides (d/dx):

= d/dx ∫( f(t) dt, 0, x) - a

= d/dx ∫( f(t) dt, 0, x) - d/dx a

Let F(t) be the anti-derivative of f(t):

= d/dx (F(x) - F(0)) - 0

= d/dx F(x) -  d/dx F(0)

F(0) is a constant.

= f(x)


Newton's Method to find the roots of f(x) can be found by the iteration:

x_(n+1) = x_n - f(x_n) / f'(x_n)

Applying that to find the roots of ∫( f(t) dt, 0, x) - a:

x_(n+1) = x_n - (∫( f(t) dt, 0, x_n) - a) / f(x_n) 


HP Prime Program INTEGRALSOLVE

Enter f(X) as a string as it will be stored in Function App variable F0.  Use X as the independent variable. 

EXPORT INTEGRALSOLVE(f,a,x)
BEGIN
// f(X) as a string, area, guess
// ∫(f(X) dX,0,x) = a
// EWS 2019-07-26
// uses Function app
LOCAL x1,x2,s,i,w;
F0:=f;
s:=0;
x1:=x;
WHILE s==0 DO
i:=AREA(F0,0,x1)-a;
w:=F0(x1);
x2:=x1-i/w;
IF ABS(x1-x2)<1 font="" then="">
s:=1;
ELSE
x1:=x2;
END;
END;

RETURN approx(x2);
END;

TI NSpire CX CAS Program INTEGRALSOLVE
(Caution:  This program needs to be typed in)

Use t as the independent variable.

Define LibPub integralsolve(f,a,x)=
Func
:© f(x), area, guess: ∫(f(t) dt,0,x = a)
:Local x1,x2,s
:s:=0
:x1:=x
:While s=0
:  x2:=x1-((∫(f,t,0,x1)-a)/(f|t=x1))
:  If abs(x2-x1)≤1E−12 Then
:    s:=1
:  Else
:    x1:=x2
:  EndIf
:EndWhile
:Return approx(x2)
:EndFunc

Examples

Example 1: 

∫( 2*t^3 dt, 0, x) = 16
Guess = 2
Root ≈ 2.3784

Example 2:

∫( sin^2 t dt, 0, x) = 1.4897
Guess = 1
(Radians Mode)
Root ≈ 2.4999

Source:

Green, Larry.  "The Second Fundamental Theorem of Calculus"  Differential Calculus for Engineering and other Hard Sciences.  Lake Tahoe Community College. http://www.ltcconline.net/greenl/courses/105/Antiderivatives/SECFUND.HTM  Retrieved July 25, 2019

Happy Solving!

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Integrals with the Error (erf) Function

Integrals with the Error (erf) Function

Introduction - Integrals Involving the Erf Function

erf(x) = 2 / √π * ∫( e^-(t^2), t, 0, x)   ** graphing calculator syntax


Problem 1:  ∫( e^-(t^2), t, a, x) with a > 0 and x > a

∫( e^-(t^2), t, 0, x) = ∫( e^-(t^2), t, 0, a)  + ∫( e^-(t^2), t, a, x)

 ∫( e^-(t^2), t, 0, x) - ∫( e^-(t^2), t, 0, a)  = ∫( e^-(t^2), t, a, x)

√π/2 * 2/√π * (  ∫( e^-(t^2), t, 0, x) - ∫( e^-(t^2), t, 0, a) ) =   √π/2 * 2/√π * ∫( e^-(t^2), t, a, x)

√π/2 * ( 2/√π *  ∫( e^-(t^2), t, 0, x) - 2/√π *  ∫( e^-(t^2), t, 0, a) ) =   ∫( e^-(t^2), t, a, x)

√π/2 * ( erf(x) - erf(a)  ) =   ∫( e^-(t^2), t, a, x)


 ∫( e^-(t^2), t, a, x) = √π/2 * (erf(x) - erf(a))

Problem 2:  ∫( e^-((t + b)^2), t, 0, x)

∫( e^-((t + b)^2), t, 0, x)

Substitutions:
u = t + b
u - b = t
du = dt

Limits:
t  = x, u = x + b
t = 0, u = b

∫( e^-((t + b)^2), t, 0, x)

= ∫( e^-(u^2), u, b, b+x)

Note √π/2  * 2 /√π  = 1

=  √π/2  * 2 /√π  * ∫( e^-(u^2), u, b, b+x)

Using problem 1 as an aid to help us get:

= √π/2  * (eft(x + b) - erf(x))

∫( e^-((t + b)^2), t, 0, x) = √π/2  * (eft(x + b) - erf(x))

Problem 3:   ∫( e^(-a*t^2 + b), t, 0, x)

∫( e^(-a*t^2 + b), t, 0, x)

= ∫( e^(-a*t^2) * e^b, t, 0, x)

e^b is a constant.

= e^b * ∫( e^(-a*t^2), t, 0, x)

Substitutions:
u = √a * t
u^2 = a * t^2

du = dt * √a
du / √a = dt

Limits:
t  = x, u = √a * x
t = 0, u = 0

= e^b * 1/√a  * ∫( e^(-u^2), u, 0, √a * x)

= e^b * 1/√a  * √π/2  * 2 /√π  * ∫( e^(-u^2), u, 0, √a * x)

= e^b * 1/√a  * √π/2  * erf(√a * x)

∫( e^(-a*t^2 + b), t, 0, x) = e^b * 1/√a  * √π/2  * erf(√a * x)

Problem 4:  ∫( e^-(ln^2 t)/t, t, 1, x)

∫( e^-(ln^2 t)/t, t, 1, x)

Substitutions:
u^2 = ln^2 t
u = ln t
du = 1/t dt
t  du = dt

1/t * t = 1

Limits:
 t = x, u = ln x
t = 1, u = ln 1 = 0

= ∫( e^-(u^2), u, 0, ln x)

= √π/2  * 2 /√π  * ∫( e^-(u^2), u, 0, ln x)

= √π/2  * erf(ln x)


∫( e^-(ln^2 t)/t, t, 1, x) = √π/2  * erf(ln x)

Problem 5:  2/√π * ∫(e^(t^2), t, 0, x)

2/√π * ∫(e^(t^2), t, 0, x)

Substitutions:
-u^2 = t^2
u^2 = -t^2
u = i * t   where i = √ -1
du = i dt
- i du = dt  because 1/i = -i

Limits:
t = x, u = i *x
t = 0, u = 0

= 2/√π * ∫(e^-(u^2) * -i, t, 0, i*x)

= 2/√π * -i * ∫(e^-(u^2), t, 0, i*x)

= -i * erf(i*x)

2/√π * ∫(e^(t^2), t, 0, x) = -i * erf(i * x)

Alternative substitution:
-u^2 = t^2
i * u = t
u = -i * t
du = -i dt
i du = dt

Limits:
t = x, u = -i *x
t = 0, u = 0

Then we get

2/√π * ∫(e^(t^2), t, 0, x)

= 2/√π * i * ∫(e^-(u^2), t, 0,- i*x)

= i * erf(-i * x)

In summary....

2/√π * ∫(e^(t^2), t, 0, x) = -i * erf(i * x) = i * erf(-i * x)

These are several of integrals using the erf function. - Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

HP 41/DM41L and TI-60X: Exponentiation of Large Numbers

HP 41/DM41L  and TI-60X:  Exponentiation of Large Numbers

But Why a Program when we have Button?

This is true.  What this program does is allow for calculation of y^x when results in answers greater than 9.999999999 * 10^9.  The number is broken up into the form:

mantissa * 10^exponent

Let n = y^x.  Then:

n = y^x

Taking the logarithm of both sides:

log n = log (y^x)
log n = x log y

A number can be split into its fractional and integer part:

log n = frac(x log y) + int(x log y)

Take the antilog of both sides:

n = 10^( frac(x log y) + int(x log y) )
n = 10^( frac(x log y) ) * 10^( int(x log y) )

where
mantissa = 10^( frac(x log y) )
exponent = int(x log y)

HP 41/DM 41L Program BIGPOW

Input:
Y stack:  y
X stack:  x

Output:
Y:  mantissa (shown first)
X:  exponent

01 LBL T^BIGPOW
02 X<>Y
03 LOG
04 *
05 ENTER↑
06 FRC
07 10↑X
08 STOP
09 X<>Y
10 INT
11 RTN

TI-60 Program:  Big Powers

Input:
Store y in R1 and x in R2

Output:
R1 = mantissa (shown first), R2 = exponent

(Step,  Key Number, Key)
00, 71, RCL
01, 02, 2
02, 65, *
03, 71, RCL
04, 01, 1
05, 43, log
06, 95, =
07, 61, STO
08, 02, 2
09, 78, Frac
10, 12, INV
11, 43, log
12, 61, STO 
13, 01, 1
14, 13, R/S
15, 71, RCL
16, 02, 2
17, 79, Intg
18, 13, R/S
19, 22, RST

Examples

Example 1:  25^76.   y = 25, x = 76

Result: 
Mantissa = 1.75162308
Exponent = 106
25^76 ≈ 1.75162308 * 10^106

Example 2:  78^55.25,  y = 78, x = 55.25

Result:
Mantissa = 3.453240284
Exponent = 104
78^55.25 ≈ 3.543240284 * 10^104

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

HP Prime and HP 41C/DM 41L: Sum of Two Squares

HP Prime and HP 41C/DM 41L:  Sum of Two Squares 

Introduction

Given a positive integer n, can we find two non-negative integers x and y such that:

n = x^2 + y^2

(x and y can be 0, n is assumed to be greater than 0)

There are several theorems and lemmas that are connected to this famous problem.  As a point of interest, I will briefly describe them here. 

1.  n does not have a representation (n can't be written as x^2 + y^2) if any of n's prime factors is congruent to 3 mod 4 and is raised to an odd power.

2.  If n has a representation, then for an integer k, k^2*n also has a representation.

3.  If n is prime and congruent to 1 mod 4, then n has a representation.  (n has the form of n = 4w + 1 for some non-negative integer w).

The program presented here is the use of iterations to find all possible pairs which fit n = x^2 + y^2.   Some integers do not have representations, others have more than one.  The program will show all possible combinations. 

HP Prime Program SUM2SQ

EXPORT SUM2SQ(n)
BEGIN
// EWS 2019-07-21
// breaking n into a sum of 2 squares
LOCAL r,j,k,l;
// we can more than 1 representation
r:=IP((n/2)^0.5);
l:={};
FOR j FROM 0 TO r DO
k:=(n-j^2)^0.5;
IF FP(k)==0 THEN
l:=CONCAT(l,
{STRING(j)+"^2 + "+
STRING(k)+"^2 = "+
STRING(n)});
END;
END;

RETURN l;
END;


HP 41C/DM 41L Program SUMSQRS

Registers  used:
R00 = n
R01 = counter
R02 = temporary

01 LBL T^SUMSQRS
02 FIX 0
03 STO 00
04  2
05  /
06  SQRT
07  INT
08  1000
09  /
10  STO 01
11  LBL 00
12  RCL 00
13  RCL 01
14  INT
15  X↑2
16  -
17  SQRT
18  STO 02
19  FRC
20  X=0?
21  GTO 01
22 GTO 02
23 LBL 01
24 RCL 01
25 INT
26 T^X = 
27 ARCL X
28  AVIEW
29  STOP
30  RCL 02
31  T^Y = 
32 ARCL X
33 AVIEW
34 STOP
35 LBL 02
36  ISG 01
37  GTO 00
38  T^END
39  VIEW
40  FIX 4
41  RTN

Examples

Example 1:  n = 325
325 = 1^2 + 18^2
325 = 6^2 + 17^2
325 = 10^2 + 15^2

Example 2:  n = 530
530 = 1^2 + 23^2
530 = 13^2 + 19^2

Source:

Dudley, Underwood.  "Elementary Number Theory"  2nd Ed. Dover Publications: New York.  1978. ISBN 978-0-486-46931-7

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

TI Nspire CX II: Distributing Money Among Friends

TI Nspire CX II:  Distributing Money Among Friends

The Problem

This post is inspired from a Instagram post by @gercekboss at www.eylemmath.weebly.com.  The problem states that:

$100 is to be distributed among 100 friends:
1% is given to the 1st friend
2% is given to the 2nd friend
3% is given to the 3rd friend
and so on.

I have a TI Nspire CX II Document that covers this problem, which can be downloaded here:  https://drive.google.com/file/d/1aewZ-8wU2JjZxgGyBpNHivEYLU5QAt31/view?usp=sharing

Name of the document:  distribute percent.tns 

Using the TI Nspire CX II to set up two sequences to model this problem as two sequences:

Balance Sequence (Problem 1.3, plot is in green)

u1(n) = u1(n-1) - u1(n-1) * n%
Intial Terms: 100
0 ≤ n ≤ 99 step 1

Amount Sequence (Problem Page 1.7, plot is in red)

u2(n) = u1(n-1) *n%
Initial Terms: 0
0 ≤ n ≤ 99 step 1

By analysis, the 10th friend will get the biggest distribution under this plan, $6.28.  By the time we get to the 38th friends, the distribution becomes meaningless because the calculated amount will be less than one penny.

Regression Analysis

I performed curve fitting analysis for the balance and amount distributed on problems page 1.5 and 1.9, respectively.

Balance can be estimated by the logistic equation:

balance(n) ≈ 112.84/(1 + 0.1*e^(0.23n)) - 0.54  for 1 ≤ n ≤ 40

The maximum absolute error is 0.67. 

Amount can be estimated by the quartic polynomial:

amount(n) ≈ -4.01E-5*n^4 + 4.31E-3 * n^3 - 0.15 * n^2 + 1.89 * n - 1.26
for 1 ≤ n ≤ 40

I did slightly better, the maximum absolute error is 0.52.

Generalizing the Problem

Now, let's take any money of money and any amount of friends for whom to distribute that money.  The rate can increase per friend by at different rates:

a*n% + b%

If you want 2%, 4%, 6%, etc, let a = 2, b = 0

If you want each friend to get 1/2% more than the last, let a = 1/2, b = 0.

If you want the first friend to get 2%, second to get 3%, etc, let a = 1, b =  1.

The sequences are set up as:

Balance Sequence

u1(n) = u1(n-1) - u1(n-1) * (a*n% + b%)
Intial Terms: c
0 ≤ n ≤ l step 1

Amount Sequence

u2(n) = u1(n-1) *(a*n% + b%)
Initial Terms: 0
0 ≤ n ≤ l step 1

Variables:
a = rate parameter
b = rate parameter
c = initial amount of money
l = number of friends

You can customize the problem in Problem 2 of the document.  Change the variables on Problem Page 2.2 and see the results on Problem Pages 2.3 and 2.4. 

Source:

Question 251. eleymmath.    https://eylemmath.weebly.com/algebra/100-question  Retrieved July 21, 2019. 

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

TI-84 Plus CE: Rectangular Waveguide - Cutoff Frequency

TI-84 Plus CE:  Rectangular Waveguide - Cutoff Frequency 

Introduction

The program WAVEGDE computes the cutoff frequency of a rectangular guide.  The waveguide has energy propagating in the TE_m,n mode with dimensions a x b (in meters) with relative permittivity E.  The cutoff frequencies is calculated as:

V = speed of light / √E
where speed of light = 299,792,458 m/s

Cutoff frequency (in Hz):
f_c = V/2 * √( (m/a)^2 + (n/b)^2 )

The second part of the calculation will depend on the target frequency, either the attenuation or phase and group velocities will be calculated:

If f < f_c, attenuation (in nerpers) is calculated:

attenuation = (2 π f_c)/V * √( 1 - (f/f_c)^2 )

If f > f_c, the phase and group velocities are calculated:

phase velocity = V / √( 1 - (f_c/f)^2 )

group velocity = V * √(1 - (f_c/f)^2 )

TI-84 Plus CE Program WAVEGDE

"EWS 2019-07-14"
Disp "RECTANGULAR WAVEGUIDE"
Input "WIDE (M): ",A
Input "NARROW (M): ",B
Input "PERMITTIVITY: ",E
Disp "HALF-WAVE VARIATIONS", "T-M,N"
Input "M: ",M
Input "N: ",N
Input "FREQUENCY (HZ):",F
299792458/√(E)→V
V/2*√((M/A)²+(N/B)²)→C
Disp "CUTOFF FREQ: (HZ)",C
Pause 
If F
Then
(2πC)/V*√(1-(F/C)²)→U
Disp "ATTENUATION:",U,"NEPERS"
Else
V/√(1-(C/F)²)→P
V*√(1-(C/F)²)→G
Disp "PHASE VELOCITY: M/S",P
Disp "GROUP VELOCITY: M/S",G
End

Examples

Example 1:
A rectangular waveguide has the dimensions 2.7 cm x 1.09 cm with permittivity of 1.    (a = 2.7/100, b = 1.09/100, E = 1).  The waves propagate in the TE_2,1 mode  (m = 2, n = 1).  The target frequency is 10 GHZ (10E9 Hz).

Results:
Cutoff Frequency = 1.767490017E10 Hz
Attenuation = 305.4488993 nepers

Example 2: 
Same as Example 1, but set target frequency as 20 GHZ.

Results:
Cutoff Frequency = 1.767490017E10 Hz
Phase Velocity = 640,624,938.6 m/s
Group Velocity = 140,293,504.8 m/s

Source:

"Rectangular Waveguide Calculations"  HP-65 E.E. Pac 2  Hewlett Packard  Cupertino, CA  1977

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Solving a System of Conic Equations: An Unorthodox Method

Solving a System of Conic Equations: An Unorthodox Method 

Introduction

Problem:  Find the solutions to the following systems of equations:

ax^2 + by^2 = c
rx^2 + sy^2 = t

Assume that:
Only a or b can be negative, while c > 0, and
Only r or s can be negative, while t > 0.

Can use the matrix method to solve for x and y?

[ [a, b], [r, s] ] * [ [ x^2 ], [ y^2 ] ] = [ [ c ], [ t ] ]

[ [ x^2 ], [ y^2 ] ] = [ [a, b], [r, s] ]^(-1) * [ [ c ], [ t ] ]

If x^2 and y^2 are both positive, the conic curves have intersecting points.

Example 1 - Two Ellipses:

3x^2 + 5y^2 = 10
4x^2 + y^2 = 6

roots: 
x^2 = 1.17647058824,  x = ± 1.0846522891
y^2 = 1.29411764706,  y = ± 1.1375291799

Intersection Points:
( 1.0846522891, 1.1375291799 ), ( -1.0846522891, 1.1375291799 ),
( 1.0846522891, -1.1375291799 ), (- 1.0846522891, -1.1375291799 )

Example 2 - An Ellipse and a Hyperbola:

3x^2 + 5y^2 = 10
4x^2 -  y^2 = 6

roots:
x^2 = 1.73913043478, x = ± 1.31876094679
y^2 = 0.956521739126, y = ± 0.978019293849

Intersection Points:
( 1.31876094679, 0.978019293849 ), ( 1.31876094679, -0.978019293849 ),
( -1.31876094679, 0.978019293849 ), ( -1.31876094679, -0.978019293849 )

Example 3 - An Ellipse and a Hyperbola, part II:

3x^2 + 5y^2 = 10
-4x^2 + y^2 = 6

roots:
x^2 = -.86956217393
y^2 = 2.51273913043

There are no intersection points.

The program SYSCONIC (HP Prime) solves the system.  Output:  a 6 x 2 matrix:

[ [ x,   0 ]
[  y,  0 ]
[  a*x^2 + b*y^2,  c ]
[  r*x^2 + s*y^2, t ]
[  a*(-x)^2 + b*(-y)^2,  c ]
[  r*(-x)^2 + s*(-y)^2, t ] ]

The last four rows are to check solutions.

EXPORT SYSCONIC()
BEGIN
HComplex:=11; // complex
// 2019-07-15 EWS
// System conic equations
LOCAL a,b,c,r,s,t;
LOCAL m0,m1,m2,m3,x,y;
INPUT({a,b,c,r,s,t},
"ax^2+by^2=c  rx^2+sy^2=t",
{"a: ","b: ","c: ",
"r: ","s: ","t: "});
m0:=[[a,b],[r,s]];
m1:=[[c],[t]];
m2:=m0^(−1)*m1;
x:=m2(1,1);
x:=√x;
y:=m2(2,1);
y:=√y;
MSGBOX(x);
MSGBOX(y);
// results and check
m3:=MAKEMAT(0,6,2);
m3(1,1):=x;
m3(2,1):=y;
m3(3,1):=a*x^2+b*y^2;
m3(3,2):=c;
m3(4,1):=r*x^2+s*y^2;
m3(4,2):=t;
m3(5,1):=a*(−x)^2+b*(−y)^2;
m3(5,2):=c;
m3(6,1):=r*(−x)^2+s*(-y)^2;
m3(6,2):=t;
END;

Remarks

It is a wise idea to check the answers to make sure that they are correct, this is kind of an unorthodox approach. 

If we tried the matrix method on a system of equations like:

x^2 + 4x - y = 2
-3x^2 + x - 2y = 5

the method will not work.  Each term must have a different variable. 


Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

A Choice of Discounts

A Choice of Discounts

Situation

We are ready to purchase two items.  Let A be the cost of the more expensive item, while B is the cost of the less expensive item (A > B).  We have two coupons:

Coupon 1:  Take 10% off of all purchases  (10% off the cost of both A and B)

Coupon 2:  Take 20% off of the most expensive item  (20% off of A, B remains at full price)

The store allows for only one coupon to be used per transaction.  Which coupon gives us the best benefit?  Let's assume that we are only buying two items.

We will not worry about sales tax in this situation.

Calculating the Costs After Discounts

Using Coupon 1:  Let T1 be the total cost using coupon 1.

T1 = (100% -10%) * (A + B)
T1 = 90% * (A + B)
T1 = 0.9 * (A + B)

Using Coupon 2:  Let T2 be the total cost using coupon 2.

Here, A is more expensive, so:

T2 = (100% -20%) * A + B
T2 = 80% * A + B
T2 = 0.8 * A + B

In summary, the total costs are:

T1 = 0.9 * (A + B)
T2 = 0.8 *A + B

Comparing the Two Coupons

Is there a price point of both A and B where the total costs of the two items is the same?  This occurs when T1 = T2.

T1 = T2

0.9 * (A + B) = 0.8 * A + B
0.9 * A + 0.9 * B = 0.8 * A + B
0.1 * A + 0.9 * B = B
0.1 * A = 0.1 * B
A = B

The costs are the same when the cost of both items are the same.

But we assume that A > B.    So, either T1 will be more or T2 will be more.  To determine this, we will need to perform a subtraction.  Let's pick T1 - T2.

If T1 -T2 > 0, then T1 costs more and using the second coupon, 20% of the more expensive item,  is preferable.

If T1 - T2 < 0, then T2 costs more and using the first coupon, 10% of both items, is preferrable.

T1 - T2
= 0.9 * A + 0.9 * B - (0.8 * A +  B)
= 0.1 * A - 0.1 * B
= 0.1 * (A - B)

Since A represents the cost of the more expensive item, A > B, and A - B > 0.  Therefore 0.1 * (A - B) > 0 and T1 > 0.   The second coupon is preferable.

Case Studies

A = $29.99,  B = $24.99:  T1 = $49.48,  T2 = $48.98  (not much difference)

A = $29.99,  B = $13.99:  T1 = $39.58,  T2 = $37.98

A = $29.99, B = $5.99;  T1 = $32.38,  T2 = $29.98  (difference between T1 and T2 increase as discrepancy of costs increase)

Hopefully this helps on all of our future shopping trips,

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

TI-84 Plus CE and HP Prime: Calculating Color Temperature

TI-84 Plus CE and HP Prime:  Calculating Color Temperature

Introduction

We can estimate a color temperature for a given color.  However, the estimate is most useful when the color is close to white or the black body curve.  The color temperature, known as the CCT, is just one particular aspect of a color and its light quality. 

For more detailed information, please check the article by Waveform Lighting, https://www.waveformlighting.com/tech/calculate-color-temperature-cct-from-cie-1931-xy-coordinates.  The article will be listed in the Sources section at the end of this blog entry.

Calculation

The program RGBTEMP will estimate CCT from a color given its RGB (red-green-blue) coordinates. 

The calculation involves several steps:

1.  Scaling RGB values from a scale of 0-255 to 0-1.  This is accomplished by dividing each of the parameters of the color's RGB value by 255.  Standard RGB values (sRGB) are used. 

2.  The scaled RGB values will be converted to CIE 1931 XYZ color space values.  In 19031, the International Commission on Illumination created the CIE 1931 color space, designed to describe how the distribution of color wavelengths affect how colors are perceived.  The conversion is done in two parts.

2a.  Adjust each of the scaled RGB values to account for the gamma correction:

Let I be the scaled parameter for R, G, and B.  Then:

If I ≤ 0.04045, Then:

I_adj = I/12.92

Else:

I_adj = ( (I + 0.055) / 1.055) ^ 2.4

2b.  Calculate the XYZ values.   This can be accomplished by matrix multiplication:

[ [ X ] [ Y ] [ Z ] ] =

[[ 0.4124, 0.3576, 0.1805 ] [ 0.2126, 0.7152, 0.0722 ] [ 0.0193, 0.1192, 0.9504 ]]
*  [ [ R ] [ G ] [ B ] ]

3.  With the XYZ coordinates, calculate CIE Yxy coordinates.  Y is already done, so only (small) x and (small) y are needed:

x = X / (X + Y + Z)

y = Y / (X + Y + Z)

4.  Calculate the CCT.  There are several methods: one of them is a cubic approximation from Calvin McCamy:

n = (x - 0.3320) / (0.1858 - y)

CCT = 437 * n^3 + 3601 * n^2 + 6861 * n + 5517

CCT is in Kelvins.

TI-84 Plus Program RGBTEMP

"2019-07-11 EWS"
Disp "COLOR TEMP FROM RGB"
Input "RED   :",R
Input "GREEN :",G
Input "BLUE  :",B

[[R/255][G/255][B/255]]→[I]

For(I,1,3)
If [I](I,1)≤0.04045
Then
[I](I,1)/12.92→[I](I,1)
Else
(([I](I,1)+.055)/1.055)^2.4→[I](I,1)
End
End

[[.4124,.3576,.1805][.2126,.7152,.0722][.0193,.1192,.9504]]*[I]→[J]

[J](1,1)/([J](1,1)+[J](2,1)+[J](3,1))→X
[J](2,1)/([J](1,1)+[J](2,1)+[J](3,1))→Y

(X-.332)/(.1858-Y)→N
437*N^3+3601*N^2+6861*N+5517→C
Disp "CCT (K):",C

HP Prime Program RGBTEMP

The HP Prime version returns CCT, the XYZ coordinates, and the x and y parameters in a four element list.

EXPORT RGBTEMP(R,G,B)
BEGIN
// 2019-07-11 EWS
// RGB to color temperature
LOCAL M0,M1,X,Y,N,C;
LOCAL I,U,V,N,CCT;

// initialize
M0:=[[R/255],[G/255],[B/255]];

// correction
FOR I FROM 1 TO 3 DO
IF M0(I,1)≤0.04045 THEN
M0(I,1):=M0(I,1)/12.92;
ELSE
M0(I,1):=((M0(I,1)+0.055)
/1.055)^2.4;
END;
END;

// RGB to CIE XYZ 1931
M1:=[[0.4124,0.3576,0.1805],
[0.2126,0.7152,0.0722],
[0.0193,0.1192,0.9504]]*M0;

// XYZ to Yxy
X:=M1(1,1)/(M1(1,1)+M1(2,1)+
M1(3,1));
Y:=M1(2,1)/(M1(1,1)+M1(2,1)+
M1(3,1));

// xy to CCT (Kelvins)
// McCamy approximation
N:=(X-0.3320)/(.1858-Y);
C:=437*N^3+3601*N^2+6861*N+5517;

RETURN {C,M1,X,Y};

END;

Examples
(CCT rounded to four decimal places)

Remember the closer to the Black Body Curve, the more accurate the answer. 

Red:  RGB (255, 0, 0)
CCT:  3034.8988 K

Green:  RGB (0, 255, 0)
CCT:  6068.7576 K

Orange:  RGB (255, 165, 0)
CCT:  2429.5395 K

Sky Blue:  RGB (136, 206, 235)
CCT:  13207.1056 K

White:  RGB (255, 255, 255)
CCT:  6506.6551 K

Sources

"Convert color data into different standards and color spaces" and "Color math and programming code examples"   EASYRGB.  IRO Group Limited.  2019. https://www.easyrgb.com/en/convert.php  and https://www.easyrgb.com/en/math.php    Retrieved July 11, 2019

"Calculate color temperature (CCT) from CIE 1931 xy coordinates"  Waveform Lighting. 2019.  https://www.waveformlighting.com/tech/calculate-color-temperature-cct-from-cie-1931-xy-coordinates   Retrieved July 6, 2019

"sRGB"  Wikipedia.  Last edited May 18, 2019.  https://en.wikipedia.org/wiki/SRGB
Retrieved June 16, 2019.

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Algebra: Solving for an Immediate Value To Calculate A Desired Result

Algebra:  Solving for an Immediate Value To Calculate A Desired Result

Introduction

Considering the following problem:

We have to equations:

F = f(X)
A = a(X)

Assume that X is close to zero (0) as possible.  Therefore, we are only considering one root.

We are given the value of F and we are tasked to find A.  In order to solve for A, we must solve for F first.

Example 1:

f(X) = 2 *X
a(X) = 1/(X^3 + 3)
F = -0.04

Solving for X:
-0.04 = 2 * X
-0.02 = X

Calculating A:
A = 1/((-0.02)^3 + 3)
A = 0.3333342222

Example 2:

f(X) = e^(X^2 - 0.5)
a(X) = ln(X + 1)/(ln X + 1)
F = 10

10 = e^(X^2 - 0.5)
ln 10 = X^2 -0.5
X^2 = ln 10 + 0.5
(we'll only use the principal root in this example)
X = √(ln 10 + 0.5)
X ≈ 1.674092319

Substituting X into a(X) to get:
A ≈ 0.6491313602

See you next time,

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

TI-84 Plus CE: Solving Algebraic Equations Step by Step

TI-84 Plus CE:  Solving Algebraic Equation Scripts

Introduction

There are four programs that uses the wait and string commands to provide a step by step instruction on how to solve the following four equations:

1.  a * x + b = c
2.  a * x^n + b = c
3.  a^x = b
4,  (a + x)^2 = b

We will need the TI-84 Plus CE.  The original TI-84 Plus with the monochrome screen will not have the necessary commands.

ALGDEMO1 

"EWS 2019-05-30"
Disp "SOLVE AX + B = C"
Prompt A,B,C
(C-B)/A→X
ClrHome
Disp toString(A)+"X + "+toString(B)+" = "+toString(C)
Wait 1.5
Disp toString(A)+"X + "+toString(B)+" - "+toString(B)+" = "+toString(C)+" - "+toString(B)
Wait 1.5
Disp toString(A)+"X = "+toString(C-B)
Wait 1.5
Disp toString(A)+" X / "+toString(A)+" = "+toString(C-B)+" / "+toString(A)
Wait 1.5
Disp "X = "+toString(X)

ALGDEMO2

"EWS 2019-05-30"
Disp "SOLVE AX^N+B=C"
Prompt A,B,C,N
((C-B)/A)^(1/N)→X
ClrHome
Disp toString(A)+"X^"+toString(N)+"+"+toString(B)+"="+toString(C)
Wait 1.5
Disp toString(A)+"X^"+toString(N)+"+"+toString(B)+"-"+toString(B)+"="+toString(C)+"-"+toString(B)
Wait 1.5
Disp toString(A)+"X^"+toString(N)+"="+toString(C-B)
Wait 1.5
Disp toString(A)+"X^"+toString(N)+"/"+toString(A)+"="+toString(C-B)+"/"+toString(A)
Wait 1.5
Disp "X^"+toString(N)+"="+toString((C-B)/A)
Wait 1.5
Disp "X"+"="+toString((C-B)/A)+"^(1/"+toString(N)+")"
Wait 1.5
Disp "PRINCIPAL ROOT:","X="+toString(X)

ALGDEMO3

"EWS 2019-05-30"
Disp "SOLVE A^X=B"
Prompt A,B
log(B)/log(A)→X
ClrHome
Disp toString(A)+"^X="+toString(B)
Wait 1.5
Disp "log("+toString(A)+")^X=log("+toString(B)+")"
Wait 1.5
Disp "X log("+toString(A)+")=log("+toString(B)+")"
Wait 1.5
Disp "DIVIDE BY log("+toString(A)+")"
Wait 1.5
Disp "X = "+toString(X)

ALGDEMO4

"EWS 2019-05-30"
Disp "SOLVE (A+X)^2=B"
Prompt A,B
­√(B)-A→X
√(B)-A→Y
ClrHome
Disp "("+toString(A)+"+X)^2="+toString(B)
Wait 1.5
Disp toString(A)+"+X=+-√("+toString(B)+")"
Wait 1.5
Disp toString(A)+"+X-"+toString(A)+"=+-√("+toString(B)+")-"+toString(A)
Wait 1.5
Disp "X=+-√("+toString(B)+")-"+toString(A)
Wait 1.5
Disp "TWO ROOTS:"
Disp "X = "+toString(X)
Disp "X = "+toString(Y)

Eddie

All original content copyright, © 2011-2019.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author