How to choose a value from the plot?

 plot(contact_geometry.y-d,contact_geometry.r_L)

Could someone please tell me, is it possible to do coding in such a way that one random value of Y is chosen for each iteration automatically.

For eg, when I run the program for the first time, R should be 0.4494, next time 0.4657. And this should rely on the above plot command. Please help.

Answer 1

[Star Strider]

Not certain what you want, so here are two possibilities:

For any y on the domain (-0.03, 0.01):

f = @(x)  1./(1-75*x);          % Substitute your function
x = 0.04*rand-0.03;             % Define domain
y = f(x);                       % Calculate ‘y’

For your particular three values:

v = [0.4494  0.4542  0.4657];
y = v(randi(3));                % Randomly choose one

[Priya]

Thank you. I get this error.

Attempted to access f(-0.0131295); index must be a positive integer or logical.

[Star Strider]

When I use my function and:

f(-0.0131295)

I get:

ans =

     503.8513e-003

If your function is not also named f, you will get that error. You have to substitute your function for mine, and if you call it something other than f, you have to change the y line to call your function correctly.

My code is simply an example, since you did not post your function. You have to adapt my code to fit your situation.

[Priya]

Yes. I'm getting it now. Thank you. But I have a doubt, I get a value of y which is not in the range shown in the plot above (which is 0.445 to 0.475).

Logically if I calculate y, I should be getting a value which falls within this limit right? But I'm getting y values like 0.8580, 0.5227 etc. could you please tell me why is that so?

[Star Strider]

I will have to see your function. I made up a function that loosely approximates yours simply to demonstrate the code. I have no idea what your data are or what the function is that generated them.

Post your function, explain the variables and parameters, and tell me what you are doing. What does the function represent? Do you need to check for bounds on x or y?

[Star Strider]

This polynomial approximation gives a good fit to your data:

load('wheel_rail_AW.mat')
S = contact_geometry
d = 0.01;

x = linspace(min(S.y-d), max(S.y-d),size(S.y-d,2));
[p,ps,smu] = polyfit(S.y-d, S.r_L, 8);      % Fit polynomial
v = polyval(p, x, ps, smu);                 % Calculate values of fit for a given s

figure(1)
plot(S.y-d, S.r_L, '.-')
hold on
plot(x, v, '-r', 'LineWidth', 1.5)
hold off
grid

I kept the plot block in so you can see the effectiveness of the fit.

I saved the [p,ps,smu] information in ‘contact_geom_poly.mat’ that I’ve attached here. All you need to do is load it, then provide an x value for this line:

y = polyval(p, x, ps, smu);

and get the y value. If you give it any value of x in the interval (-0.03,0.01), it should give you reasonably accurate values for y.

There are other approaches for this of course, one being interpolation, but that is significantly more complicated.

EDIT —

This is the easiest way to use the polynomial I fit to your data:

load('contact_geom_poly.mat')  
f = @(x) polyval(p, x, ps, smu);

y1 = f(-0.018)                      % Scalar argument
y2 = f([-0.0021  0.0067])           % Vector argument

yields:

y1 =

     449.4902e-003

y2 =

     454.3047e-003   465.7238e-003

These compare almost exactly to the data on the plot in your Question.

[Priya]

Thanks very much for your code. I'll try, work on this a bit and will let you know.

[Star Strider]

My pleasure!

[Priya]

It's all working now, but I expect 'y' to generate a random value for each 'x' instead of these three points alone shown in the plot above. I have just shown it for an example.

So each time I run the program, y should automatically show a different value within the mentioned x and y-axis limits.

Thanks

[Star Strider]

I seem not to be making myself clear.

The full code I intend you to implement is:

load('contact_geom_poly.mat')  
f = @(x) polyval(p, x, ps, smu);
x = 0.04*rand-0.03;             % Define domain
y = f(x);

That’s simply part of the example code I gave you before, and the specific polynomial fit to your data that I derived.

This code should do exactly what you want!

[Priya]

Thank you so much again. I got it!

[Star Strider]

My pleasure!

I actually demonstrated that in my comment headed: ‘Star Strider on 24 May 2014 at 13:02’ about 7-8 comments above this one.

Here is some revised code to make it as easy as possible for you to calculate your regression coefficients and structures for various dependent variables. The routine does its best to fit your data, then creates ‘.mat’ files for the polynomial fit output so you can easily use them in your application:

datamatfile = 'wheel_rail_AW.mat';
dirmats = which(datamatfile);           % Find directory for input data file
dirsep = strfind(dirmats, '\');         % Get path to ‘datamatfile’
pathname = dirmats(1:max(dirsep));      % Create path (‘fileparts’ returns an empty string for some reason)
load(datamatfile)                       % Load data
S = contact_geometry;                   % Create data structure
d = 0.01;                               % Define ‘d’
SNames = fieldnames(S)                  % Variable names in ‘datamatfile’
% ———————— Prompt for dependent variable for polynomial regression ————————
vblnr = listdlg('PromptString', {'Select a dependent variable' 'for polynomial regression:'}, 'ListString', SNames, 'SelectionMode','single')
vbl = char(SNames(vblnr));              % Variable to be regressed

x = linspace(min(S.y-d), max(S.y-d),size(S.y-d,2));
[p,ps,smu] = polyfit(S.y-d, eval(['S.' vbl]), 51);      % Fit polynomial
v = polyval(p, x, ps, smu);              % Calculate values of fit for a given s

figure(1)
plot(S.y-d, eval(['S.' vbl]), '.-')
hold on
plot(x, v, '-r', 'LineWidth', 1.5)
hold off
xlabel('y-d')
ylabel('Amplitude')
title(sprintf('Plot of ‘%s’ as a function of [y-d]', vbl))
grid

% ———————————————— Save regression information and notify ————————————————
save([pathname 'contact_geom_poly_' vbl '.mat'], 'p', 'ps', 'smu')
h1 = msgbox({'Polynomial regression structures saved as' ['contact_geom_poly_' vbl '.mat']});

Your data have some significant discontinuities, so to be as robust as possible with this, I am using an excessively high-degree polynomial. It fits most of your data reasonably well, but ‘rings’ at the discontinuities. There is no way to avoid that, because ‘ringing’ will occur at the ends if it fits the discontinuities without ringing. This way, it ‘rings’ a bit at the discontinuities and at the ends, but minimally at both. This is about the best you can expect from it.

The routine searches for the directory ‘wheel_rail_AW.mat’ is in and stores the ‘.mat’ files it creates in the same directory.

-------------------------------------------------------------------------------------------------------------------------------------

NOTE — I will hear quite justifiable screams in my direction at my using a 51-degree polynomial. The purpose here is to provide a reasonable approximation of data with some saturation discontinuities over a narrow domain to make the interpolation as easy as possible in the intended application. Reading the entire thread provides the necessary context.

[Priya]

Sorry, I can see that you have clearly mentioned everything I need to do. But since I haven't done this before, I have no idea how to assign values for the regression coefficients.

Also, is the whole chunk of code above demonstrates how to create a .mat file for storing the regression coefficients? Please do this favour.

[Star Strider]

Not a problem. I just wanted to be sure you know that I’m not leaving you stranded.

You don’t need to assign the values for the regression coefficients. The polyfit function does that for you. All you have to do is choose the variable you want polyfit to provide the polynomial fit for. The code I posted in my previous comment will do everything for you.

Run the entire code in my previous comment. (You may want to make a separate script file for it.) You will see that it offers you a menu of variables to regress against (y-d), does the regression, plots it, stores the estimated parameters and structures in a ‘.mat’ file that incorporates the name of the variable at the end of the ‘.mat’ file name, and tells you what it did.

All you have to do then is run these lines in your code:

load('contact_geom_poly_r_L.mat')  
f = @(x) polyval(p, x, ps, smu);
x = 0.04*rand-0.03;             % Define domain
y = f(x);

to use the fit for ‘r_L’. Substitute the ‘.mat’ file name of your choice in the load statement for the other variables.

It will do its best to produce the y values you want for the ‘.mat’ file you choose.

[Priya]

Thanks a ton and above all thanks for your patience. It was of great help.

[Star Strider]

My pleasure!

[Priya]

Hi,

I would like to know if there's any function other than polyfit to perform the above.

[Star Strider]

You may want to experiment with the spline function. Experiment with the interp1 as well (with spline or piecewise interpolation options), and there could be some File Exchange functions, but I’m not aware of any others.

The problem with some of your data are the saturation discontinuities. You could test for them, use polyfit on the smooth portions of your data, then have your function return a constant for x values above the saturation discontinuity. I used polyfit because it was the logical first choice for what you were originally doing, and gave decent fits, even at the discontinuities.

[Star Strider]

Nothing actually wrong. It’s just that in that region, sin(x) = x. Change the line computing the sine to:

 y=sin(100*x);

and you’ll see that it works as intended.

Change the plot command to:

 plot(x,y,'.',xi,yi, 'or')

so that the one interpolated point is more visible.

[Star Strider]

I have no idea. I’m lost. What did you change?

[Star Strider]

It would seem to me that:

Y = f_fun(X);

is essentially equivalent to (and can be replaced by):

yi=interp1(x,y,xi);

You don’t need the polynomial fit if you’re using interp1.

[Star Strider]

To force the plot to only display y-axis limits only between 0.445 and 0.475, add this to your plot statements:

axis([xlim    0.445  0.475])

[Star Strider]

It works. There are simply not many data in that region:

load('wheel_rail_AW.mat')
S = contact_geometry
x = linspace(min(S.y), max(S.y),size(S.y,2));
y=sin(100*x);
xi = 0.04*rand-0.03;
yi=interp1(x,y,xi);
figure(1)
plot(x,y,'.',xi,yi,'or')
axis([xlim    0.445  0.475])

[Star Strider]

It seems you can.

The ‘NaN’ is most likely caused by a value for xi that is outside of the x data range. A solution for that is to change your interp1 call to:

yi = interp1(x_new, y_new, xi, 'linear', 'extrap');

That won’t change its behaviour (the default method is 'linear'), except to allow linear extrapolation.

[Priya]

Now yi happens to be out of range

[Star Strider]

That is to be expected. You are extrapolating.

If you want to choose random numbers only on your ‘x_new’ range:

xi = min(x_new) + (max(x_new)-min(x_new)) * rand;

[Priya]

Perfect! Thanks once again.

[Star Strider]

Again, my pleasure!