h = @(x1,x2,x3) x1.*(t - x2).^x3;
g = @(x1,x2,x3) exp(-0.5.*((y - h(x1,x2,x3))./sigma).^2 );
f = @(x1,x2,x3,theta) g(x1,x2,x3).*theta(1).*theta(2).*theta(3).*theta(4);
J1 = @(theta) integral3(@(x1,x2,x3) f(x1,x2,x3,theta), 0, 1, 0, t, 0.3, 1);
When you use integral3(), the function must accept three arrays the same size, and must return an output the same size. It must perform element-wise multiplication.
So integral3 is going to pass arrays for x1, x2, x3, and @(x1,x2,x3) f(x1,x2,x3,theta) must return something the same size. The result of integral3() will be a scalar.
This has the direct implication that the size of the output of @(x1,x2,x3) f(x1,x2,x3,theta) must not depend upon the size of theta. And it also means that
J2 = @(theta) sum(log(J1(theta)));
the J1 is going to return a scalar, and it does not make sense to take sum() of a scalar.
Your code is obviously hoping that you can use integral3() to return something the same size as theta, but integral3() cannot do that.
To return something the same size as theta, you are going to have to use arrayfun, such as
J2 = @(Theta) sum(log(arrayfun(J1, Theta)));
Or perhaps it is not theta you are concerned about...
Look again at the sequence. If y were non-scalar then g would not return something the same size as the x values, and that would mean that f would return something that was not the same size as the x values, and that would mean that @(x1,x2,x3) f(x1,x2,x3,theta) would return something that was not the same size as the x values, but returning something the same size is a requirement of integral3() . So y cannot be a non-scalar.
Perhaps you were hoping that integral3() would return one value for each y... but it cannot do that. You would have to arrayfun() over the Y values and pass the current y into the functions.
