Showing posts with label inverse. Show all posts
Showing posts with label inverse. Show all posts

Sunday, September 3, 2023

Python: Basic Functions with Quaternions

 Python:   Basic Functions with Quaternions   



Calculating with Quaternions



Quaternions are four-dimensional hypercomplex numbers that is written as the form:


r = a * i + b * j + c * k + d


We can also write quaternions as a four dimensional vector [a, b, c, d].


In computer graphics, quaternions are used to draw, display, move, and rotate 3D objects.  The constant term, d, is a mass-point, while the other three terms describe the location of vector point.



Some Mathematical Properties of Quaternions



The conjugate of a quaternion is similar to what a conjugate of a regular complex number would look like:


r* =  - a * i - b * j - c * k + d  


Similarly, the norm of a quaternion works in a familiar fashion:


| r | = abs(r) = √(a^2 + b^2 + c^2 + d^2)


We can multiply a quaternion by any scalar, s, as follows:


s * r = s*  a * i + s * b * j + a* r3 * k + a* r4  


To find the unit vector (unit quaternion), multiply the quaternion by the scalar 1/|r|.


Adding quaternions is communicative and associative.  However, multiplying quaternions is associative, but not communicative.  In short, when it comes to multiplying quaternions, order matters.  


For the complex coefficients:


i^2 = j^2 = k^2 = -1


i * j =  k,    j * i = -k 


j * k = i,   k * j = -i


k * i = j,  i * k = -j




List of Functions Included in QUAT.py



Quaternions are entered as vectors:


r = a * i + b * j + c * k + d -> [a, b, c, d]


qplus(r, q):   addition, r + q

qminus(r, q):  subtraction, r - q

qtimes(r,q):  multiplication,  r * q

qdivide(r,q):   division (multiply by the reciprocal):  r / q  = r * q^-1

qnorm(r):  norm or magnitude

qconj(r):  conjugate

qscale(r, a):   multiply quaternion r by a scalar a

qinv(r):  inverse or reciprocal of a quaternion,  1/r = r* / |r|

qsquare(r):  r^2


qhelp():  brings up the list of functions in QUAT plus.  


To use these functions on other Python scripts:


from QUAT import *


Important QUAT should also import the math module.  




Python:  Quaternions


This should work on all Python platforms since only the math module is used.  This was programmed on a TI-84 Plus CE Python calculator.




Python Code:  QUAT.py


# arithmetic with quaternions

# ews 2023-06-20


from math import *


# help

def qhelp():

  s="Quaternions - qhelp()"

  # add \n for new lines

  s+="\nEnter r=ai+bj+ck+d as [a,b,c,d]."

  s+="\nqplus(r,q) +  qminus(r,q) -"

  s+="\nqtimes(r,q) *  qdivide(r,q) /"

  s+="\nqnorm(r) magnitude \nqconj(r) conjugate \nqscalar(r,a) scalar mult. by a"

  s+="\nqinv(r) r**-1 qsquare(r) r**2"

  # need print to look nice

  return print(s)


# add

def qplus(r,q):

  return [r[i]+q[i] for i in range(4)]


# subtract

def qminus(r,q):

  q=[-i for i in q]

  return qplus(r,q)


# conjugate

def qconj(r):

  return [-r[i] for i in range(3)]+[r[3]]


# norm

def qnorm(r):

  s=sum([i**2 for i in r])

  return sqrt(s)


# inverse/reciprocal

def qinv(r):

  r=qconj(r)

  s=qnorm(r)

  s*=s

  return [i/s for i in r]


# scalar multiplication

def qscalar(r,a):

  return [a*i for i in r]


# multiplication

def qtimes(r,q):

  s=[0,0,0,0]

  s[0]=r[0]*q[3]+r[1]*q[2]-r[2]*q[1]+r[3]*q[0]

  s[1]=-r[0]*q[2]+r[1]*q[3]+r[2]*q[0]+r[3]*q[1]

  s[2]=r[0]*q[1]-r[1]*q[0]+r[2]*q[3]+r[3]*q[2]

  s[3]=r[3]*q[3]-r[0]*q[0]-r[1]*q[1]-r[2]*q[2]

  return s


# division

def qdivide(r,q):

  q=qinv(q)

  return qtimes(r,q)


# qsquare

def qsquare(r):

  return qtimes(r,r)


qhelp()



Download the script here:  https://drive.google.com/file/d/1iNDxw82Gm1nOVJ7jDHdnIIxA4tcaV0p3/view?usp=sharing



Source


Ron Goldman,

Understanding quaternions,

Graphical Models,

Volume 73, Issue 2,

2011,

Pages 21-49,

ISSN 1524-0703,

https://doi.org/10.1016/j.gmod.2010.10.004.

(https://www.sciencedirect.com/science/article/pii/S1524070310000172)


Science Direct requires a paid subscription or access through a university.


Alternate source:  https://www.researchgate.net/publication/220632454_Understanding_quaternions




All original content copyright, © 2011-2023.  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. 


Thursday, September 8, 2022

Casio fx-991EX Classwiz Tips: Normal Distribution

Casio fx-991EX Classwiz Tips:  Normal Distribution


This week I am going to show some things that can be done with the Casio fx-991EX Classwiz.  


The Normal Distribution 


The Classwiz's Distribution mode has calculations for the following probability distributions:


1.  Normal, f(x) = exp(-1/2 * ((x - μ) / σ)^2 ) / (σ * √(2 * π))

2.  Binomial

3.  Poisson


This is Mode 7.


Normal CD - Finding the Area


[ MENU ], 7: Distribution, 2: Normal CD


The CDF function calculates the area (probability) between two limits.   The lower and upper tail areas require the limits to be -∞ and +∞, respectively.  However, the Classwiz does not provide values for -∞ and +∞.  For the best estimate, I suggest using -7 and 7 are the limits.  The calculator calculates the probability from -7 to 7 to be 1.


Example:

(For all the examples, the standard values μ = 0, σ = 1)


Lower tail area to x = 3:  lower limit = -7, upper limit = 3; area:  0.9986501019


Lower limit = -2, upper limit = 2; area:  0.954499736


Upper tail from x = -1:  lower limit = 1, upper limit = 7; area:  0.1586552539


Inverse CD


[ MENU ], 7: Distribution, 3: Inverse Normal


This calculation gets the point value for a lower tail distribution (-∞ to x).


Example:

(For all the examples, the standard values μ = 0, σ = 1)


Area = 0.5; xInv:  0


Area = 0.6; xInv: 0.2533470931


Area = 0.7; xInv: 0.5244004382


Eddie


All original content copyright, © 2011-2022.  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. 


Saturday, November 21, 2020

Numworks: 3 x 3 Matrices

 Numworks:   3 x 3 Matrices


The script invthree.py calculates the inverse and determinant of a 3 x 3 matrix.  


Matrices:


[[ a1,  a2,  a3 ]

[ b1, b2, b3 ]

[ c1, c2, c3 ]]



Numworks script invthree.py:

(701 bytes)


from math import *


# 2020-11-12 EWS

print("3x3 Matrix Inverse")

print("[[a1,a2,a3]")

print("[b1,b2,b3]")

print("[c1,c2,c3]]")

a1=float(input('a1: '))

a2=float(input('a2: '))

a3=float(input('a3: '))

b1=float(input('b1: '))

b2=float(input('b2: '))

b3=float(input('b3: '))

c1=float(input('c1: '))

c2=float(input('c2: '))

c3=float(input('c3: '))

# determinant

d=a1*(b2*c3-b3*c2)-a2*(b1*c3-b3*c1)+a3*(b1*c2-b2*c1)

# inverse

d1=(b2*c3-c2*b3)/d

d2=-(a2*c3-c2*a3)/d

d3=(a2*b3-a3*b2)/d

e1=-(b1*c3-c1*b3)/d

e2=(a1*c3-c1*a3)/d

e3=-(a1*b3-a3*b1)/d

f1=(b1*c2-b2*c1)/d

f2=-(a1*c2-a2*c1)/d

f3=(a1*b2-a2*b1)/d

print("det=")

print(d)

print("inv=")

print([d1,d2,d3])

print([e1,e2,e3])

print([f1,f2,f3])


Numworks page:  https://workshop.numworks.com/python/ews31415/invthree


Example:


[[ -5.4, 3.3, -1.7 ], [ 0.6, 8.3, 5.3 ] [ 5.5, 5.4, 1.9 ]] 


returns 


[ -0.05493...,  -0.06604..., 0.13508... ]

[ 0.11974... ,  -0.00389..., 0.11798... ]

[ -0.18130..., 0.20224..., -0.20006... ]


Source:


wikiHow Staff "How to Find the Inverse of a 3x3 Matrix"  WikiHow.   Last Updated November 5, 2020.  https://www.wikihow.com/Find-the-Inverse-of-a-3x3-Matrix  Retrieved November 12, 2020.


Eddie


All original content copyright, © 2011-2020.  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. 


Thursday, September 12, 2019

HP 32SII and TI-66: Curve Fitting

HP 32SII and TI-66:  Curve Fitting

Introduction

The curve fitting program uses the linear regression module to determine the parameters b ("intercept") and m ("slope") in non-linear curves using following transformations:

Logarithmic Regression:  y = b + m * ln x
Transformations:  ( ln x, y, b, m )

Inverse Regression:  y = b + m / x
Transformations:  ( 1/x, y, b, m )

Exponential Regression:  y = b * e^(m * x)
Transformation:  ( x, ln y, e^b, m )

Power Regression:  y = b * x^m
Transformation:  ( ln x, ln y, e^b, m )

Geometric (Exponent) Regression:  y = b * m^x
Transformation:  ( x, ln y, e^b, e^m )

Simple Logistic Regression:  y = 1 / (b + m * e^(-x))
Transformation:  ( e^(-x), 1/y, b, m )

HP 32SII Program:  Curve Fitting

Note:
1.  This can be adapted into the HP 35S under one label.  Just take note of the where the label points are.
2.  The total amount of bytes used is 90.
3.  Flags 1 and 2 are used.  If flag 1 is set, e^m is calculated as slope.  If flag 2 is set, e^b is calculated as intercept.

Program:
// Initialize - LBL X
LBL X
CF 1
CF 2
CLΣ
0
RTN

// Calculation - LBL Y
LBL Y

FS? 2
e^x
STO B
VIEW B
m
FS? 1
e^x
STO M
VIEW M

STO R
VIEW R
RTN

// Logarithmic Regression - LBL L
LBL L
LN 
R/S
GTO L

// Inverse Regression - LBL I
LBL I
1/x
R/S
GTO I

// Exponential Regression - LBL E
LBL E
SF 2
x<>y
LN 
x<>y
R/S
GTO E

// Power Regression - LBL P
LBL P
SF 2
LN 
x<>y
LN 
x<>y
R/S
GTO P

// Geometric/Exponent Regression - LBL G
LBL G
SF 1
SF 2
x<>y
LN
x<>y
R/S 
GTO G

// Simple Logistic Regression - LBL S
LBL S
+/-
e^x
x<>y
1/x 
x<>y
STOP 
GTO S

Instructions:
1.  Clear the statistics data and flags by pressing [XEQ] X.
2.  Enter data points, run the proper label, and press [ Σ+ ] or [ Σ- ].

For example, for Logarithmic fit:
y_data [ENTER] x_data [XEQ] L [ Σ+ ]

Subsequent Data:
y_data [ENTER] x_data [R/S] [ Σ+ ]

This scheme allows for undoing data:
y_data [ENTER] x_data [XEQ] L [ Σ- ]

3.  Calculate intercept (B), slope (M), and correlation (R), press [XEQ] Y.

TI-66 Program:  Curve Fitting

Notes:
1.  This program should be able to entered on a TI-58, TI-58C, or TI-59.  At the time of the posting, I have not done it, so I don't have the key codes.
2.  94 steps are used.  [INV] [SBR] is merged into the RTN step.
3.  Flags 1 and 2 are used.  If flag 1 is set, e^m is calculated as slope.  If flag 2 is set, e^b is calculated as intercept.

Program:
// Initialize - key [ A ]
000 LBL
001 A
002 INV
003 ST.F
004 01
005 INV
006 ST.F
007 02
008 CSR
009 0
010 RTN

// Calculation - key [ A' ]
011 LBL
012 A'
013 OP
014 12
015 INV
016 IF.F
017 02
018 (   // left parenthesis
019 INV
020 LN X
021 LBL
022 (  // left parenthesis
023 STO
024 08
025 R/S
026 X<>T
027 INV
028 IF.F
029 01
030 )  // right parenthesis
031 INV
032 LN X
033 LBL
034 )  // right parenthesis
035 STO 
036 07
037 R/S
038 OP
039 13
040 STO
041 09
042 RTN

// Logarithmic Regression - key [ B ]
043 LBL
044 B
045 LN X
046 X<>T
047 R/S
048 RTN

// Inverse Regression - key [ C ]
049 LBL 
050 C
051 1/X
052 X<>T
053 R/S
054 RTN

// Exponential Regression - [ D ]
055 LBL 
056 D
057 ST.F
058 02
059 X<>T
060 R/S
061 LN X
062 R/S
063 RTN

// Power Regression - [ B' ]
064 LBL
065 B'
066 ST.F
067 02
068 LN X
069 X<>T
070 R/S
071 LN X
072 R/S
073 RTN

// Geometric/Exponent Regression - [ C' ]
074 LBL
075 C'
076 ST.F
077 01
078 ST.F
079 02
080 X<>T
081 R/S
082 LN X
083 R/S
084 RTN

// Simple Logistic Regression - [ D' ]
085 LBL
086 D'
087 +/-
088 INV
089 LN X
090 X<>T
091 R/S
092 1/X
093 R/S
094 RTN

Instructions:
1.  Clear the statistics data and flags by pressing [  ].
2.  Enter data points: enter x, run the proper label, enter y, press [R/S] and press [2nd] ( Σ+ ) or [INV] [2nd] ( Σ+ )  (for  Σ- ).

For example, for Logarithmic fit:
x_data [ B ] y_data [R/S]  [2nd] (Σ+)

This scheme allows for undoing data:
x_data [B] y_data [R/S] [INV] [2nd] (Σ+)

3.  Calculate intercept (B), slope (M), and correlation (R), press [2nd] [ A' ].

Examples

All results are rounded.

Example 1: Logarithmic Regression
Data (x,y):
(33.8, 102.4)
(34.6, 103.8)
(36.1, 105.1)
(37.8, 106.9)

Results:
B:  -33.4580
M:  38.6498
R:  0.9941

y ≈ -33.4580 + 38.6498 ln x

Example 2:  Inverse Regression
Data (x,y):
(100, 425)
(105, 429)
(110, 444)
(115, 480)

B:  823.80396
M:  -40664.72143
R:  -0.91195

y ≈ 823.80396 - 40664.72143/x

Example 3: Simple Logistic Regression
Data (x,y):
(1, 11)
(1.3, 9.615)
(1.6, 8.75)
(1.9, 8.158)
(2.6, 7.308)

B: 0.14675
M: -0.15487
R:  -0.99733

y ≈ 1 / (0.14675 - 0.15487*e^(-x))


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.


Saturday, November 4, 2017

HP Prime: Best Regression Fit

HP Prime:  Best Regression Fit

The program BESTFIT compares a set of regressions to determine a best fit.  This simulates a feature presented on the Hewlett Packard HP 48S, HP 48G, HP 49G, and HP 50g.  BESTFIT compares the correlations of the following four regression models:

1.  Linear:  y = a * x + b
2.  Logarithmic:  y = a * ln x + b
3.  Exponential:  y = b * a^x   (y = b * e^(ln a * x))
4.  Power: y = b * x^a

The output is a two element list:  a string of the best fit equation and its corresponding correlation.

HP Prime Program:  BESTFIT

EXPORT BESTFIT(L1,L2)
BEGIN
// 2017-11-02 EWS
// Simulate Best Fit
// HP 48GX,49G,50g

// initialize
LOCAL clist,m,c,v,s,n;
clist:={0,0,0,0};

// test
// correlation is linear only
// requires approx()
// linear y=a*x+b
clist[1]:=approx(correlation(L1,L2));
// log y=a*LN(x)+b
c:=approx(correlation(LN(L1),L2));
IF IM(c)==0 THEN
clist[2]:=c;
END;
// exponential y=b*e^(a*x)
c:=approx(correlation(L1,LN(L2)));
IF IM(c)==0 THEN
clist[3]:=c;
END;

// power y=b*a^x
c:=approx(correlation(LN(L1),LN(L2)));
IF IM(c)==0 THEN
clist[4]:=c;
END;

// test
// POS(L0^2,MAX(L0^2))
m:=POS(clist^2,MAX(clist^2));
c:=clist[m];

IF m==1 THEN
v:=linear_regression(L1,L2);
s:=STRING(v[2])+"+"+STRING(v[1])+
"*X";
RETURN {s,c};
END;

IF m==2 THEN
v:=logarithmic_regression(L1,L2);
s:=STRING(v[2])+"+"+STRING(v[1])+
"*LN(X)";
RETURN {s,c};
END;

IF m==3 THEN
v:=exponential_regression(L1,L2);
s:=STRING(v[2])+"*"+STRING(v[1])+
"^X";
RETURN {s,c};
END;

IF m==4 THEN
v:=power_regression(L1,L2);
s:=STRING(v[2])+"*X^"+
STRING(v[1]);
RETURN {v,s};
END;

END;

Examples – BESTFIT:

Example 1:

X
y
1.05
10.45
2.28
11.33
4.20
16.38
6.34
28.87

Result: {“7.78250037344*1.21713132288^X”, 0.981261724397}
Example 2:

X
Y
-10
82
-5
41
5
-42
10
-79

Result:  {“0.5+ -8.1*X”, -0.999801918343}


Example 3:

X
Y
10
2.278
11
2.666
12
2.931
13
3.212

Result:  {“-5.79464870365+3.51429873838*LN(X)”, 0.998295735284}

BESTFIT2:  An Extended Version

Version 2 adds the following regressions: 

5.  Inverse:  y = b + a/x
6.  Simple Logistic:  y = 1/(b + a*e^(-x))
7.  Simple Quadratic:  y = b + a*x^2
8.  Square Root:  y = √(a*x + b)

HP Prime Program:  BESTFIT2

EXPORT BESTFIT2(L1,L2)
BEGIN
// 2017-11-02 EWS
// Simulate Best Fit
// HP 48GX,49G,50g
// additional models

// initialize
LOCAL clist,m,c,v,s;
clist:={0,0,0,0,0,0,0,0};

// test
// correlation is linear only
// requires approx()

// linear y=a*x+b
clist[1]:=approx(correlation(L1,L2));

// log y=a*LN(x)+b
c:=approx(correlation(LN(L1),L2));
IF IM(c)==0 THEN
clist[2]:=c;
END;

// exponential y=b*e^(a*x)
c:=approx(correlation(L1,LN(L2)));
IF IM(c)==0 THEN
clist[3]:=c;
END;

// power y=b*a^x
c:=approx(correlation(LN(L1),LN(L2)));
IF IM(c)==0 THEN
clist[4]:=c;
END;

// inverse y=b+a/x
IF POS(L1,0)==0 THEN
clist[5]:=approx(correlation(1/L1,
L2));
END;

// simple logistic
IF POS(L2,0)==0 THEN
clist[6]:=approx(correlation(e^(−L1),
1/L2));
END;

// simple quadratic
clist[7]:=approx(correlation(L1^2,
L2));

// square root
IF ΣLIST(L2≥0)==SIZE(L2) THEN
clist[8]:=approx(correlation(L1,
L2^2));
END;



// test
// POS(L0^2,MAX(L0^2))
m:=POS(clist^2,MAX(clist^2));
c:=clist[m];

IF m==1 THEN
v:=linear_regression(L1,L2);
s:=STRING(v[2])+"+"+STRING(v[1])+
"*X";
END;

IF m==2 THEN
v:=logarithmic_regression(L1,L2);
s:=STRING(v[2])+"+"+STRING(v[1])+
"*LN(X)";
END;

IF m==3 THEN
v:=exponential_regression(L1,L2);
s:=STRING(v[2])+"*"+STRING(v[1])+
"^X";
END;

IF m==4 THEN
v:=power_regression(L1,L2);
s:=STRING(v[2])+"*X^"+
STRING(v[1]);
END;

// inverse
IF m==5 THEN
v:=linear_regression(1/L1,L2);
s:=STRING(v[2])+"+"+STRING(v[1])+
"/X";
END;

// simple logistic
IF m==6 THEN
v:=linear_regression(e^(−L1),
1/L2);
s:="1/("+STRING(v[2])+"+"+
STRING(v[1])+"*e^(−X))";
END;

// simple quadratic
IF m==7 THEN
v:=linear_regression(L1^2,L2);
s:=STRING(v[2])+"+"+STRING(v[1])+
"*X^2";
END;

// square root
IF m==8 THEN
v:=linear_regression(L1,L2^2);
s:="√("+STRING(v[2])+"+"+
STRING(v[1])+"*X)";
END;

RETURN {s,c};
END;

Examples – BESTFIT2:

Example 4:

X
y
1.05
10.45
2.28
11.33
4.20
16.38
6.34
28.87

Result:  {“9.0573970192+0.48023219108*X^2”, 0.994491728382}

Example 5:

X
y
1
8.4853
2
8.9443
4
9.7980
7
10.9546

Result: {“√(63.9995089183+8.00048914762*X”, 0.999999999759}

Example 6:

X
y
1
1
2
-0.5
3
-1
4
-1.25

Result:  {“-2+3/X”, 1}

Eddie

This blog is property of Edward Shore, 2017

HP 20S: Acoustics Programs

HP 20S: Acoustics Programs Program A: Speed of Sound in Dry Air cs = 20.05 × √(273.15 + T°C) Code: 01: 61, 41, A: LBL A 02: ...