Multiple regression: Difference between revisions

=={{header|Perl 6}}==

We're going to solve the example on the Wikipedia article using [https://github.com/grondilu/clifford Clifford], a [https://en.wikipedia.org/wiki/Geometric_algebra geometric algebra] module. Optimization for large vector space does not quite work yet, so it's going to take (a lof of) time and a fair amount of memory, but it should work.

Let's create four vectors containing our input data:

<math>\begin{align}

\mathbf{w} & = w^k\mathbf{e}_k\\

\mathbf{h_0} & = (h^k)^0\mathbf{e}_k\\

\mathbf{h_1} & = (h^k)^1\mathbf{e}_k\\

\mathbf{h_2} & = (h^k)^2\mathbf{e}_k

\end{align}</math>

Then what we're looking for are three scalars <math>\alpha</math>, <math>\beta</math> and <math>\gamma</math> such that:

<math>\alpha\mathbf{h0} + \beta\mathbf{h1} + \gamma\mathbf{h2} = \mathbf{w}</math>

To get for instance <math>\alpha</math> we can first make the <math>\beta</math> and <math>\gamma</math> terms disappear:

<math>\alpha\mathbf{h0}\wedge\mathbf{h1}\wedge\mathbf{h2} = \mathbf{w}\wedge\mathbf{h1}\wedge\mathbf{h2}</math>

Noting <math>I = \mathbf{h0}\wedge\mathbf{h1}\wedge\mathbf{h2}</math>, we then get:

<math>\alpha = (\mathbf{w}\wedge\mathbf{h1}\wedge\mathbf{h2})\cdot\tilde{I}/I\cdot\tilde{I}</math>

'''Note:''' a number of the formulae above are invisible to the majority of browsers, including Chrome, IE/Edge, Safari and Opera. They may (subject to the installation of necessary fronts) be visible to Firefox.

<lang perl6>use Clifford;

my @height = <1.47 1.50 1.52 1.55 1.57 1.60 1.63 1.65 1.68 1.70 1.73 1.75 1.78 1.80 1.83>;

my @weight = <52.21 53.12 54.48 55.84 57.20 58.57 59.93 61.29 63.11 64.47 66.28 68.10 69.92 72.19 74.46>;

my $w = [+] @weight Z* @e;

my $h0 = [+] @e[^@weight];

my $h1 = [+] @height Z* @e;

my $h2 = [+] (@height X** 2) Z* @e;

my $I = $h0∧$h1∧$h2;

my $I2 = ($I·$I.reversion).Real;

say "α = ", ($w∧$h1∧$h2)·$I.reversion/$I2;

say "β = ", ($w∧$h2∧$h0)·$I.reversion/$I2;

say "γ = ", ($w∧$h0∧$h1)·$I.reversion/$I2;</lang>

{{out}}

<pre>α = 128.81280357844

β = -143.1620228648

γ = 61.960325442</pre>

This computation took over an hour with the april 2016 version of rakudo on MoarVM, running in a VirtualBox linux system guest hosted by a windows laptop with a i7 intel processor.

=={{header|Phix}}==

{{trans|ERRE}}

<lang Phix>constant N = 15, M=3

sequence x = {1.47,1.50,1.52,1.55,1.57,

1.60,1.63,1.65,1.68,1.70,

1.73,1.75,1.78,1.80,1.83},

y = {52.21,53.12,54.48,55.84,57.20,

58.57,59.93,61.29,63.11,64.47,

66.28,68.10,69.92,72.19,74.46},

s = repeat(0,N),

t = repeat(0,N),

a = repeat(repeat(0,M+1),M)

for k=1 to 2*M do

for i=1 to N do

s[k] += power(x[i],k-1)

if k<=M then t[k] += y[i]*power(x[i],k-1) end if

end for

end for

-- build linear system

for row=1 to M do

for col=1 to M do

a[row,col] = s[row+col-1]

end for

a[row,M+1] = t[row]

end for

puts(1,"Linear system coefficents:\n")

pp(a,{pp_Nest,1,pp_IntFmt,"%7.1f",pp_FltFmt,"%7.1f"})

for j=1 to M do

integer i = j

while a[i,j]=0 do i += 1 end while

if i=M+1 then

?"SINGULAR MATRIX !"

?9/0

end if

for k=1 to M+1 do

{a[j,k],a[i,k]} = {a[i,k],a[j,k]}

end for

atom Y = 1/a[j,j]

a[j] = sq_mul(a[j],Y)

for i=1 to M do

if i<>j then

Y=-a[i,j]

for k=1 to M+1 do

a[i,k] += Y*a[j,k]

end for

end if

end for

end for

puts(1,"Solutions:\n")

?columnize(a,M+1)[1]</lang>

{{out}}

<pre>

Linear system coefficents:

{{ 15.0, 24.8, 41.1, 931.2},

{ 24.8, 41.1, 68.4, 1548.2},

{ 41.1, 68.4, 114.3, 2585.5}}

Solutions:

{128.8128036,-143.1620229,61.96032544}

</pre>