Iterative Proportional Fitting

[Stefan Hrafn Jonsson]

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
-------------------------------------------------------------