On Sun, 22 Jan 2012, Barton Willis wrote:
> Your change of variable converts the lower limit of integration to
> either -infinity for infinity (depending on the sign of B). The
> function changevar doesn't handle such things correctly, I think.
> You could try changing the lower limit to a new symbol---I wasn't
> able to push this any further.
> For integer N, one thing that might work:
> (%i1) assume(C>0, B>0);
> (%o1) [C>0,B>0]
> (%i2) marginal : integrate((1/t-C)^(N-1)*exp(A-B/t^2),t,0,1/C);
> (%o2) integrate(%e^(A-B/t^2)*(1/t-C)^(N-1),t,0,1/C)
> (%i3) makelist(ev(marginal, N : k, integrate),k,1,3);
> (%o3) [(gamma_incomplete(-1/2,B*C^2)*%e^A*sqrt(B))/2,-(gamma_incomplete(-1/2,B*C^2)*%e^A*sqrt(B)*C-gamma_incomplete(0,B*C^2)*%e^A)/2,-(-gamma_incomplete(-1/2,B*C^2)*%e^A*B*C^2+2*gamma_incomplete(0,B*C^2)*%e^A*sqrt(B)*C-gamma_incomplete(1/2,B*C^2)*%e^A)/(2*sqrt(B))]
Thanks Barton. In the end, I did the binomial expansion by hand, and
the substitution piecemeal, like this:
declare(N,integer) ;
assume(N > 0) ;
assume(N > 1) ;
assume(B > 0) ;
marginalshouldbe : 'integrate((1/t-C)^(N-1)*exp(A-B/t^2),t,0,1/C) ;
theintegrand : (1/t-C)^(N-1)*exp(A-B/t^2) ;
theintegrandexpanded : (1/t)^(N-1)*sum(gamma(N)*(-C)^j/(gamma(j+1)*gamma(N-j)*t^(N-j-1)),j,0,N-1)*exp(A-B/t^2) ;
theintegrandtermj : (1/t)^(N-1)*gamma(N)*(-C)^j/(gamma(j+1)*gamma(N-j)*t^(N-j-1))*exp(A-B/t^2) ;
subst(t = B^(1/2)/u^(1/2), theintegrandtermj) ;
%*diff(B^(1/2)/u^(1/2),u) ;
integrate(%,u,inf,B*C^2) ;
marginal : sum(%,j,0,N-1) ;
--
Dan