Polynomial fit with trend1d

Q&A
1

Salut! I need help with fitting the second-degree polynomial through the data. Basically, I did this:

Blockquote
#!/bin/bash
I=“-i3,2”
J=“-JX12c”
gmt makecpt -Cjet -Ic -T2021-08-25T12/2021-09-15T12/1d > time.cpt
gmt begin wd_hv_inc png
gmt set FORMAT_DATE_MAP “jjj”
gmt basemap -R2/8/3/8 $J -Bafg -Bpx+l"Water discharge (m@+3@+s@±1@+)" -Bpy+l"Horizontal velocity (mm hr@±1@+)"
gmt plot hv_incr.dat -L -i3,1,0 -Sc5p -W1.5p,black -G+z -Ctime.cpt
gmt colorbar -DJCR -Ctime.cpt -Bsxa1D+l"Day of year"
gmt trend1d -Fxmr -i3,1 hv_incr.dat -Np2 -V > model.tmp
gmt plot model.tmp -W2p,red
gmt end show

and got multiple fitting lines:

wd_hv_inc
wd_hv_inc1888×1561 173 KB

Could you please let me know what I should tweak to get only one line of fit? I read the whole documentation on trend1d and tried many options, but failed to find a solution. The data is just a table without any discontinuities. I am using GMT 6.4.0.

I appreciate your time and help,
Anuar

2

Can you try with capital P ? -NP2
The way I understand the doc, small p gives all the possibilities up to 2.

p (polynomial with intercept and powers of x up to degree n ), P (just the single term x^n )

3

@PlanetGus, thank you very much for your answer. Unfortunately, it did not fix the problem. I feel like the problem is in a way how model.tmp file is created by trend1d. Are there any ways to control the output of this file?

4

If you get multiple lines then I suspect your input is a multiple segment file and you get one solution per segment? No other way to get more than one predicted value per x value.

5

@pwessel, thanks al lot for your response! Yes, that was my first thought but the data as I said has no discontinuities, for example with a character like this “>”. I have attached the data file and here is trend1d messages:

Blockquote
trend1d [INFORMATION]: Solving for 3 terms using a last squares norm. The model:
trend1d [INFORMATION]: y(x) = a + b * x^1 + c * x^2
trend1d [INFORMATION]: In the above, x = 2*(x - x_mid)/(x_max - x_min) for polynomials and x = 2pi(x - origin)/length for Fourier components
trend1d [INFORMATION]: The complete polynomial will be represented via Chebyshev polynomials
trend1d [INFORMATION]: Processing input table data
trend1d [INFORMATION]: Reading Data Table from file hv_incr.dat
trend1d [INFORMATION]: Read 561 data with X values from 2.908 to 7.108
trend1d [INFORMATION]: N_model Rank Chi_Squared Significance
trend1d [INFORMATION]: 3 3 0.24904146642 1
trend1d [INFORMATION]: Final model stats: N model parameters 3. Rank 3. Chi-Squared: 0.24904146642
trend1d [INFORMATION]: Model Coefficients (Polynomial): 6.49997025994 -0.775282082905 0.113014014508
trend1d [INFORMATION]: Writing Data Table to Standard Output stream

hv_incr.dat (29.9 KB)

6

This is the relevant line. A MWE should post only it plus a plot of the result

gmt trend1d -Fxmr -i3,1 hv_incr.dat -Np2 -V > model.dat

Indeed, it seems to be a bug as model points are not sorted as they should.
And when I wrote this I found the problem. Your input points are not monotonic increasing and GMT doesn’t complain.

7

Thank you @Joaquim! After sorting the x axis in ascending order before feeding to trend1d, there are no longer multiple lines. Thank you all for your time and help!!!

8

FYI, we have updated the master repo so if m is selected amount -F settings then we automatically now sort output on increasing x.

1 Like