help-octave
[Top][All Lists]

## Iterative Proportional Fitting

 From: Stefan Hrafn Jonsson Subject: Iterative Proportional Fitting Date: Thu, 26 Oct 2000 15:26:33 -0400 (EDT)

```Hi all

I am new to Octave

I am trying to find a way to use Iterative Proportional Fitting to solve
simoltanioysly 13 equations.

I do have a Maple code to do this but I would really like to do this in Octave,
partly because I have not been able to find Maple for unix on this (Penn State)
campus (except too old versions) and bacause people around tell me Octave
should  be good for this task.

So far I have been able to construct in octave all the given values needed to
solve the equations.
I am stuck in the octave code wehre I go about to solve equations. To
demostrate what equation I want I provide the Maple code used to solve this.

Msd = 0.00068
Msm = 0.13446
Mmd = 0.00037
Mmw = 0.00116
Mmv = 0.06563
Mwd = 0.00148
Mwm = 0.03334
Mvd = 0.00113
Mvm = 0.5128
Ms = [Msd+Msm,0,0,0;-Msm,Mmd+Mmw+Mmv,-Mwm,-Mvm; 0,-Mmw,Mwd+Mwm,0;
0,-Mmv,0,Mvd+Mvm ]
II = [1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1]
Msplusi = II + 0.5*Ms
Msplusinv = inverse(Msplusi)
Msminusi = II - 0.5*Ms
Spi = Msminusi * Msplusinv

# Observed rates for from 1995 US MSLT, Females 22-23
aMsd=0.00058
aMsm=0.09772
aMmd=0.0003
aMmw=0.00112
aMmv=0.04588
aMwd=0.0014
aMwm=0.05119
aMvd=0.000753
aMvm=0.22632

AM =
[aMsd+aMsm,0,0,0;-aMsm,aMmd+aMmw+aMmv,-aMwm,-aMvm;0,-aMmw,aMwd+aMwm,0;0,-aMmv,0,aMvd+aMvm]

AMplusi = II + 0.5*AM

AMplusinv =inverse(AMplusi)
AMminusi  = II - 0.5*AM
Api  = AMminusi*AMplusinv

x = [76411;20716;64;1814]
y = [69250;27035;85;2432]

Cpi = Api * x

Sctot1 = Spi(1,1)+Spi(2,1)+Spi(3,1)+Spi(4,1)
Sctot2 = Spi(1,2)+Spi(2,2)+Spi(3,2)+Spi(4,2)
Sctot3 = Spi(1,3)+Spi(2,3)+Spi(3,3)+Spi(4,3)
Sctot4 = Spi(1,4)+Spi(2,4)+Spi(3,4)+Spi(4,4)

# Now how do we solve the 13 equations

#function Z = f (x)
#  Z = [z11,0,0,0;z21,z22,z23,z24;z31,z32,z33,z34;z41,z42,z43,z44]
#endfunction
#

#The Maple way of diong this.

Z :=
matrix(4,4,[z11,0,0,0i],[z21,z22,z23,z24],[z31,z32,z33,z34],[z41,z42,z43,z44]]);

solve ({ y[1,1] = z11*x[1,1],
y[2,1] = z21*x[1,1]+z22*z[2,1]+z23*x[3,1]+z24*x[4,1],
y[3,1] = z31*x[1,1]+z32*z[2,1]+z33*x[3,1]+z34*x[4,1],
y[4,1] = z41*x[1,1]+z42*z[2,1]+z43*x[3,1]+z44*x[4,1],
Sctot2/Sctot1 = (z32+z32+z42)/(11+z21+z31+z41),
Sctot3/Sctot1 = (z33+z33+z43)/(11+z21+z31+z41),
Sctot4/Sctot1 = (z34+z34+z44)/(11+z21+z31+z41),
z21*z42/(z22*z41) = Spi[2,1]*Spi[4,2]/(Spi[2,2]*Spi[4,1]),
z21*z43/(z23*z41) = Spi[2,1]*Spi[4,3]/(Spi[2,3]*Spi[4,1]),
z21*z44/(z24*z41) = Spi[2,1]*Spi[4,4]/(Spi[2,4]*Spi[4,1]),
z31*z42/(z32*z41) = Spi[3,1]*Spi[4,2]/(Spi[2,2]*Spi[4,1]),
z31*z43/(z33*z41) = Spi[3,1]*Spi[4,3]/(Spi[2,3]*Spi[4,1]),
z31*z44/(z34*z41) = Spi[3,1]*Spi[4,4]/(Spi[2,4]*Spi[4,1])},
{z11,z21,z31,z41,z22,z32,z42,z23,z33,z43,z24,z34,z44})

Any help or hints will be apreciated.

Thanks
Stefan

_________________________

Stefan Hrafn Jonsson
701 Oswald Tower
Department of Sociology and
The Population Research Institute
The Pennsylvania State University,
University Park, PA 16802
_________________________

-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

```