Matt J's solution got me started on the right path. As it turns out, the bsxfun command was unnecessary. I've posted my solution on the file exchange for anyone else who wishes to use it. Thanks to all for your help.
polyfit along a given dimension
18 次查看(过去 30 天)
显示 更早的评论
The mean function allows the user to specify a dimension along which to calculate the mean. Can I do something similar with polyfit? I have a big 3D data set. It's essentially a stack of gridded values that have evolved through time, and I'd like to make a map of the linear trend without using a loop.
0 个评论
采纳的回答
更多回答(3 个)
Matt J
2014-4-3
编辑:Matt J
2014-4-3
Assuming you only want the fit and don't want any of the error analysis output that polyfit normally gives, you can do things like
[m,n,p]=size(data3D);
data2D=reshape(data3D,m,[]);
V=bsxfun(@power,(0:m-1).'/m, [1,0] );
fit =reshape( V*(V\data2D), m,n,p); %linear fit along dim=1
And, of course, you can permute() the data3D array to achieve the same along other dims.
4 个评论
Image Analyst
2014-4-3
You can make a visualization that is a movie where go through slice by slice. I'm not sure what your trend is. Is it the mean value of each slice? Then maybe you can take the mean of each slice with mean2() and fit it to a quadratice or whatever with polyfit() and then plot it as a 1D curve with plot(). I'm not really sure what the input to polyfit would be. What are you regressing? The slice means? Can you explain more?
3 个评论
Image Analyst
2014-4-19
What is windData? Is that the 3D array? So you have 180 rows by 360 columns by 12600 frames? OK, fine. But I still don't know what linear trend is. I'm just going to have to repeat my questions. Trend of what? Is it the trend of something as a function of frame number? If so, of what? The mean? The variance? If so, get the mean and variance of each frame into an array and then run polyfit on that.
Paul Shepherd
2014-4-18
Hi Chad! I tried to send you an e-mail to make contact as well. Hopefully you got my contact info.
So, if you want a poor-mans linear fit, you might take the diff() along the appropriate dimension, and then consider the mean and variance of the diff. Would that be useful?
I don't know how often you are sampling, but it seems like the linear fit along the time axis (12600 samples) might not be the most interesting thing to look at. What if you down-sampled this to look at the average value over a longer sample period? It's been a while since I looked at decimation filters, but this might let you get the time data down quite a bit. Unfortunately, the filter() function still assumes 1-D data, so you might need to code up an algorithm to scale things up to work with the MxN timeslice as a single "data point".
Fianlly, have you looked at any parallel processing options to maybe make the brute-force method less problematic?
Hope this helps, Paul Shepherd
2 个评论
Paul Shepherd
2014-4-18
I think you can significantly improve the processing time of a full polyfit by reordering your indices. A little experiment:
x = [1, 2; 3, 4];
y = [5, 6; 7, 8];
z = zeros( 2, 2, 2);
z(:,:,1) = x;
z(:,:,2) = y;
x(1:end)
z(1:end)
Doing this, you can see that Matlab stores things column-wise, even for an N>2 D matrix. Try re-ordering your data so that the matrix is not MxNxP, where you want to do the fit across P, but MxPxN, where the outer for loop counts N, and the inner for loop counts M. That should allow the optimal memory utilization in your program, and might lead to a significant speed up in your algorithm.
John D'Errico
2014-4-22
编辑:John D'Errico
2014-4-22
Paul - The mean of diff is a very poor mans solution, when there are better answers available.
For example, consider this simple case:
y = [0:5, -3:1];
>> polyfit(1:length(y),y,1)
ans =
-0.22727 2.2727
>> mean(diff(y))
ans =
0.1
See that the mean of the diffs gives a slope where the sign is not even the same as that given by polyfit.
Sorry, but I'd avoid this idea. Again, there are trivial solutions to make the polyfit no more difficult than a matrix multiply.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!