Cauchy distribution from two zero-mean normals

See Papoulis, pp. 196–199.
PDF of bivariate normal distribution, X and Y
In[]:=
fXY=PDF[MultinormalDistribution[{{sx^2,rsxsy},{rsxsy,sy^2}}],{x,y}]
Out[]=
1
2
-
y(rsyx-sxy)
-1+
2
r
sx
2
sy
-
x(-syx+rsxy)
-1+
2
r

2
sx
sy

2π
2
sx
2
sy
-
2
r
2
sx
2
sy
In[]:=
$Assumptions=Element[{x,y,z,sx,sy,r},Reals]&&sx>0&&sy>0&&r^2<=1
Out[]=
(x|y|z|sx|sy|r)∈&&sx>0&&sy>0&&
2
r
≤1
PDF of Z = X/Y
In[]:=
fZ=2Integrate[yfXY/.x->yz,{y,0,∞}]
Out[]=
2
1-
2
r
sxsy
2π
2
sx
-4πrsxsyz+2π
2
sy
2
z
In[]:=
fZ=Simplify[fZ]
Out[]=
1-
2
r
sxsy
π
2
sx
-2πrsxsyz+π
2
sy
2
z
This is not the form given by Papoulis, but we can check it against his, which is
In[]:=
fZpap=Sqrt[1-r^2]sxsy/Pi/(sx^2(1-r^2)+sy^2(z-rsx/sy)^2)
Out[]=
1-
2
r
sxsy
π(1-
2
r
)
2
sx
+
2
sy
2
-
rsx
sy
+z
In[]:=
Simplify[fZ==fZpap]
Out[]=
True
Now let’s see if we can set up Mathematica to give us this from the CauchyDistribution function
In[]:=
cauchy=Simplify[PDF[CauchyDistribution[z0,γ],z]]
Out[]=
γ
π(
2
z
-2zz0+
2
z0
+
2
γ
)
In[]:=
cauchy/.{γ->Sqrt[1-r^2]sx/sy,z0->rsx/sy}
Out[]=
1-
2
r
sx
πsy
2
r
2
sx
2
sy
+
(1-
2
r
)
2
sx
2
sy
-
2rsxz
sy
+
2
z
In[]:=
Simplify[%]
Out[]=
1-
2
r
sxsy
π
2
sx
-2πrsxsyz+π
2
sy
2
z
In[]:=
%==fZ
Out[]=
True
So the shape parameter is
In[]:=
Sqrt[1-r^2]sx/sy
Out[]=
1-
2
r
sx
sy
And the location parameter is
In[]:=
rsx/sy
Out[]=
rsx
sy